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Abstract 

An efficient, approximate analysis for calculating 
complete three-dimensional stress fields near regions of 
geometric discontinuities in laminated composite struc- 
tures is presented. An approximate three-dimensional 
local analysis is used to determine the detailed local re- 
sponse due to far-field stresses obtained from a global 
two-dimensional analysis. The stress results from the 
global analysis are used as traction boundary conditions 
for the local analysis. A generalized plane deformation 
assumption is made in the local analysis to reduce the so- 
lution domain to two dimensions. This assumption al- 
lows out-of-plane deformation to occur. The local 
analysis is based on the principle of minimum comple- 
mentary energy and uses statically admissible stress 
functions that have an assumed through-the-thickness 
distribution. Examples are presented to illustrate the ac- 
curacy and computational efficiency of the local analy- 
sis. Comparisons of the results of the present local 
analysis with the corresponding results obtained from a 
finite element analysis and from an elasticity solution are 
presented. These results indicate that the present local 
analysis predicts the stress field accurately. Computer 
execution-times are also presented. The demonstrated 
accuracy and computational efficiency of the analysis 
make it well suited for parametric and design studies. 

Introduction 

Advanced composite materials are being used more 
extensively for aircraft primary structural components. 
To exploit fully the benefits of composite materials for 
use in these comonents, a thorough understanding of the 
corresponding response and failure characteristics of 
these structures is required. Practical forms of these pri- 
mary structural components have many local detail fea- 
tures that cause gradients in the stresses and 
displacements near these detail features. In addition, 
some structural features have thickness discontinuities 
that cause local three-dimensional interlaminar stress 
gradients. Thickness discontinuities commonly occur at 
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flange terminations in a stiffened panel, at dropped plies 
or skin thickness changes, at the edges of fail-safe straps, 
and at the edges of lap joints. Local interlaminar stress 
gradients have been identified as being significant con- 
tributors to the failure of composite structures for all 
loading conditions. An example of such a failure occurs 
in stiffened panels that are loaded into the postbuckling 
load range. 1 * Failure can occur for these panels when the 
high local stress gradients in the interface between the 
skin and the stiffener flange have stresses that exceed the 
low interlaminar strength of the composite materials. 
These stress gradients can cause the stiffeners to separate 
from the skin. 

With few exceptions, previous analyses of lami- 
nated plate structures with thickness discontinuities have 
been conducted using displacement-based finite element 
analyses. The large number of degrees of freedom re- 
quired to obtain accurate transverse stress predictions 
with this analysis method makes the use of this method 
computationally expensive or even infeasible for struc- 
tures of practical size and with practical structural design 
details. This observation is particularly true if the analy- 
sis is included in a design or structural optimization code 
where analyses are performed repetitively or, if nonlinear 
deformations associated with postbuckling or internal 
pressure are considered. 

The present paper presents an alternative approxi- 
mate analysis method for determining the three-dimen- 
sional interlaminar stresses at a local thickness 
discontinuity in a laminated composite structure. The 
analysis is limited to linear elastic material behavior and 
geometrically linear deformations. The approximate 
analysis method is based on stress functions and the prin- 
ciple of minimum complementary energy. The method is 
similar to an approximate analysis method presented by 
Kassapoglou and DiNicola, 5,6 but uses the more refined 
stress function assumptions for each layer of a laminate 
that were used by Yin 7 for analyzing the classic free edge 
problem. In addition, the present analysis includes both 
the thick and the thin regions of a panel with a thickness 
discontinuity as shown in Fig. 1 . Only the thick region of 
the panel was modeled in the previous approximate anal- 
yses. 5,6 Including both the thick and thin panel regions 
and the refined stress assumptions for each layer provides 
a more accurate stress field prediction, and a better repre- 
sentation of the load transfer between the two regions. 
Furthermore, the present analysis can model general lam- 
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inates and is not restricted to orthotropic layers as are 
currently available elasticity solutions. 2,3,4 

The flange termination region at a skin-stiffener in- 
terface is emphasized in the present paper; however, the 
method can be applied or extended to other thickness dis- 
continuity configurations that occur in practical compos- 
ite structures. The computational efficiency of the 
analysis is demonstrated by solution execution-time re- 
sults, and the accuracy of the analysis is demonstrated by 
comparisons of results from the present analysis with re- 
sults from finite element and elasticity analyses. The re- 
sults presented herein suggest that the present 
approximate analysis will provide structural designers 
with a computationally efficient means for assessing the 
influence of structural, geometric and material parame- 
ters on local three-dimensional stress fields and for de- 
signing laminated composite structures with critical 
local detail features. 

Problem Definition and Solution Approach 

The physical problem considered has a single dis- 
continuity, or change in thickness, as shown in Fig. 1 . A 
global-local analysis approach is used for the present 
analysis to make the analysis general and applicable to a 
variety of structural configurations. The global analysis 
consists of a two-dimensional finite element analysis of 
a structure subjected to prescribed loading and displace- 
ment boundary conditions. Stress resultants obtained 
from the global analysis at a point (or points) of interest 
in the structure are applied as far-field traction boundary 
conditions to a local analysis of the thickness discontinu- 
ity region. The remainder of the present paper focuses 
on the details of the local analysis. 

It is assumed that the deformations in the local 
analysis can be represented accurately by the assumption 
of generalized plane deformations. For this class of 
problems, the loading, geometric and material proper- 
ties, and stresses and strains are independent of the lon- 
gitudinal coordinate x (see Fig. 1 for the definition of the 
coordinate system). The solution domain can therefore 
be reduced to any y -z cross-sectional plane that is re- 
mote from the ends of the structure where the stresses 
and strains are influenced by the boundary conditions, 
and a state of generalized plane deformation may not ex- 
ist. 

o 

Lekhnitskii investigated the class of deformations 
for which the stress tensor is independent of the longitu- 
dinal coordinate x , and found that the most general form 
of the displacement field for this class of problems (see 
page 104 of Ref. 8) is 

u ( x, y, z) - x (Ay + Bz + C) + U (y, z) + u 
v(x,y,z) = (-7/2) Ax-txz + V(y, z) + v' (1) 

w (x, y, z) = (-1/2) Ex' + Ojc y + W (y, z) + w' 


where A and B are the negative bending curvatures in 
the x-y and x-z planes, respectively, C is the strain 
£ x at y = 0 and z = 0 and Ox is the rotation about the 
x axis. The values of A, B, C and 0 are the same for 
all layers. The functions U, V and W are unknown com- 
ponents of the displacement, and u' v' and xv' are rigid 
body displacements. 

A model of a cross-section of the structural detail 
and the applied loads used in the local analysis are shown 
in Fig. 2. The stress resultants N and A/ y shown in Fig. 
2 and the stress resultants N , M x and M are deter- 
mined from the global analysis. The local solution do- 
main is divided into two regions, A and B, as shown in 
Fig. 2a, to account for the change in geometry at the 
thickness discontinuity. Each region is assumed to be 
sufficiently long in the y -direction so that the solution 
provided by the global two-dimensional analysis is valid 
away from the discontinuity. Independent solutions are 
developed for the two regions, and inter-region continu- 
ity conditions obtained from the variational statement are 
used to connect the two solutions. Region A includes 
both the terminated laminate (shown in Fig. 1) and the 
laminate that continues into region B. Region A consists 
of N layers separated by n planar interfaces. An adhe- 
sive or bond layer may also be included between the con- 
tinuous and terminated laminates. Region B includes 
only the laminate continued from region A, and consists 
of N s layers which are separated by n s interfaces. The 
thickness of the layers are represented by 
h(k) = z (k-\) _ z {k) and are numbered consecutively 
from the top surfaces of the laminates (see Fig. 2b). 
Mathematical layers may represent individual plies that 
are reinforced by a system of parallel fibers oriented at an 
angle 0 with respect to the x axis or they may represent 
a portion of the thickness of a ply. The layers are contin- 
uous between regions. Each layer (indicated by the sub- 
script k) is regarded as a homogeneous, monoclinic, 
linear elastic material, defined by the following constitu- 
tive relation 
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where contracted notation has been used; e.g., 

Kx’ VV V °«’ % ]7 = [ ° l ’° 2 ’ ° 3 > ° 4 ’ 0 5 ’° 6] 7 < 3 > 
The corresponding notation is used for the strains. For 
generalized plane deformation, a longitudinal strain 
can be specified. Expressing the longitudinal 
stress (J^ in terms of this strain and the remaining 
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( 8 ) 


stress components gives 

0 (t) = _ S (k) cfk) y s (k) (J = 2 ; 2 , 6) (4) 

In Eq. (4) and the equations that follow, Cartesian tensor 
summation rules are assumed unless otherwise noted. 
Substituting Eq. (4) into the right-hand side of Eq. (2) 
and replacing the first row by the above expression for 
o|k) gives 
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where 

P U k) = S ij k) ~ S ii )/S lj ) (6) 
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The general expression in Eq. (8) is specialized to the 
generalized plane deformation problem by recalling that 
the stresses in all layers are independent of the axial 
coordinate jc. As a result, only variations of stresses 
that are independent of x are considered in the analysis. 
In addition, the integration in jc is carried out for a unit 
length. Carrying out the integration in jc reduces the 
volume integral in Eq. (8) to an area integral in the y - z 
plane, and the surface integral to a contour integral in 
the same plane. Making these simplifications, and 
employing the notation in Fig. 1, Eq. (8) becomes 

5n c = Z LA s u k)a ‘ k)5a j k)dA 

k= 1 () 

- J [AuSap) +Av5a 6 ^ + Aw8g 5 ( ^]<M 

A ik ) 

l [w5<r^^ + v8o 2 ^ + vv8 o^-^dz = 0 

(9) 


Local Analytical Model 


The principle of minimum complementary energy 
is used with Lekhnitskii stress functions 8 to determine 
the stresses in the local analysis. Neglecting body forces, 
the complementary energy functional for a laminated 
body consisting of N subvolumes (layers) with volume 
V (£) , is written in Cartesian coordinates as 


N 

M I 


k=\lv 


2 V 1 J 


dV- f x u.dS (7) 

J t i u v 7 


(*) 


In the above expression, aW denotes the stresses, x. 
denotes the surface tractions, 5. denotes the prescribed 
displacements, and SjW denotes the compliance matrix 
given in Eq. (2). The superscript (&) attached to the 
variables indicates that the variable is associated with 
the k th layer. The displacements w. are prescribed 
along the boundary segment S u . The stress tensor 0- 
satisfies the equilibrium equations within each layer and 
the traction continuity conditions on inter-layer bound- 
aries, and it is compatible with the prescribed tractions 
over the boundary segment S G . The functional 
JT^must be stationary with respect to the independent 
admissible variations in the stresses, and the compatibil- 
ity equations are determined from the resulting Euler- 
Lagrange variational equations. The first variation of 
the complementary energy functional is given as 


where Aw, Av, Aw are the differences between the dis- 
placements at the ends x = 0 and jc = 1 , and A ^ 
denotes the area of the k th layer in the jc = 0 plane. 
The longitudinal stress component, is determined 
from Eq. (4). 


Laver Stress Assumptions 


The structural analysis proceeds by approximating 
the stresses within the layers in the two regions shown in 
Fig. 2. The stress assumptions for the two regions are in- 
dependent, but are identical in form. Hence, the subse- 
quent discussion is developed for one region. A generic 
layer in the analysis has thickness h K ; and the upper and 
lower surfaces are defined by z ^ ~ 1 ) and z ^ , respec- 
tively, as shown in Fig. 3. A local non-dimensional co- 
ordinate 1} = (z ~ ^ ) is defined 

for the analysis. A stress field in the k th layer satisfying 
the governing equilibrium equations may be expressed in 
terms of a pair of stress functions 8 (y, T|) and 

4*(*) (y, r|) for the generalized plane deformation as- 
sumption used in the present analysis such that 


o 2 <*> = 


02 
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Approximate forms for the stress functions are ob- 
tained by assuming a solution for the stress field within a 
typical layer using separation of variables 9 ; i.e., 

oW (>-, q) = /.(*) (y) gW (q) / = / 6 

(no sum on i) 

where fW (y) are unknown functions of the in-plane 
variable y , and gW (t|) are assumed explicit polyno- 
mial shape functions of the through-the-thickness vari- 
able r\ . The through-the-thickness polynomial functions 
presented by Yin 7 and by Rose and Herakovich 10 are 
used in the present analysis. In these formulations, 
expressions for the through-the-thickness shape func- 
tions are developed by assuming that the in-plane stress 
components G^ and G^ vary linearly through the 
thickness of a layer in accordance with classical lamina- 
tion theory. Thus, from Eq. (10), the through-the-thick- 
ness shape functions for O W in a layer are defined by 
cubic polynomials in r| , and the through-the-thickness 
shape functions for y M in a layer are defined by qua- 
dratic polynomials in r| . The assumed polynomial 
stress functions for each layer (k) used in the present 
analysis are expressed as 

<t (t) (>’,q) = (1 -3q 2 + 2q 3 )F*(y) (12) 

+ (q-2q 2 + q 3 )h^G k (y) 

+ (3q 3 - 2q 3 ) F k _ , (>’) 

+ (q 3 -q 2 )/,(*>G*_,O0 


defined by z = z ( ,V) and the top surface defined by 
z = yields 

( 15 ) 

d>(D ( Z = z (°)) = 0W( Z = Z (N)) = 0 

<D,Z (1) (z = 2 ( ° 1 ) = <&, Z (Ar> (z = =0 

*KC) (z = z<°>) = ¥ (W(z = z M) =0 

In addition, free surface conditions are imposed along 
surface AB in Fig. 2 and inter-region continuity condi- 
tions are imposed along surface BC. These conditions 
are outlined in the next section. 

Application of the Variational Principle 

The unknown functions F k (y) , G^(y) , 

^(k) ( y ) * an< * H(k) (>’) ^ determined by applying Eq. 
(9) to a generic laminate with edge boundaries at 
y = Tj and y = Y 2 (see Fig. 3). In the minimization 
of the complementary energy, the second and third terms 
in Eq. (9) contribute only to the particular solutions of 
the unknown functions. Therefore, these terms are 
dropped in the formulation of the complementary solu- 
tion. Substituting for the stresses in terms of the stress 
functions in the first term in Eq. (9) and integrating by 
parts with respect to y yields 



+ (P23°,yyzz + P 3 3 % + , yyz > «<*> 


'PW(yll) = d-tl 2 )I t W (13) 

+ ( T l 2 ) L k _ i (y) + (q-q 2 )h^H k {y) 


where F k (y) , G k (y) ,L k (y) and H k (y) are unknown 
functions of the in-plane variable y . 

In addition to satisfying the equilibrium equations 
within a layer, the stresses must satisfy the traction con- 
tinuity conditions between layers, and the surface bound- 
ary conditions. Continuity of the interlaminar normal 
stress and the interlaminar shearing stresses at the k th in- 
terface z = z(k) requires that along the interface 


o(*+ O 


*<» <i>,^ +1) =o, 2 (t) 

q/(* + 1) = 4 /(*) 


(14) 


Free surface conditions are also imposed along the top 
and bottom surfaces of the laminates. For pressure- 
loaded panels, these surfaces may be exposed to normal 
pressure tractions. However, the magnitude of the pres- 
sure load should be small relative to the stresses that 
develop within the material, and it is ignored in the anal- 
ysis. Imposing the free surface conditions at the top and 
bottom surfaces of a laminate with the bottom surface 


+ x£ up 23 ®«*|),,®., y+ ivf,) 8 ®, 

k= 1 

-<l i 23®,«^J®,,„'>P»' 1 ',„> S ® 

+ <P.4®,„ *(>«'•',> 5®., 

+ (P4 5 ®.„+p 5! ®,)5® >*ii> 

= 0 

(16) 

In Eq. (16) the superscripts (k) have been dropped for 
conciseness and a comma denotes partial differentiation. 
The double integral in Eq. (16) contains the compatibil- 
ity equations associated with the variations in the stress 
tensor, and the remaining terms in the equation are the 
continuity conditions between the regions. Equation 
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(16) can be expressed in terms of the 3 n unknown func- 
tions F k (y) , G k ( y ) and L k (y) and the n + 1 
unknown functions H k (y) by substituting the stress 
function assumptions of Eqs. (12) and (13) into Eq. (16) 
and integrating through the thickness of the laminate. 
Since explicit piecewise continuous functions have been 
assumed for the through-the-thickness variations of the 
stresses, the through-the-thickness integration is per- 
formed by summing the individual layer integrals in 
z (k \ Carrying out these manipulations gives 

+ X}dy (17) 

= 0 


where the vectors {X} and {Y} have the form 
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r X, = F, 


*„ +1 = G , 


X 2n = G „ 

*2„ + . = L l 
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3 n n 
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*3„ + /V = 


{Y} = 


*1 = 


Y = F 

n n 
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y 2 n = G „ 
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^ = L n 
y 3n+l = H l 


Y 3n + N = H N 


x 3n + V + 1 
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dF 
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r 4n + /V- 

ac 


4 n+N+ 1 


f 5n + V 


ay 


ac 

n 

ay 


(18) 


Solution Approach 

The matrices [ U] , [ V] , and [W] in Eq. (17) 


have elements that are real-valued constants and are de- 
termined by evaluating the double integral term in Eq. 
(16). The matrices [P] , [Q] , [/?] , and [5] are ob- 
tained from the other terms in Eq. (16). The individual 
components of these matrices have been determined us- 
ing the symbolic manipulation computer code MACSY- 
MA. 1 1 These components are extremely lengthy and are 
not presented. Satisfying Eq. ( 1 7) for the arbitrary vari- 
ations 5 F k > 6 G k , 5 L k and 6 H k yields a system of lin- 
ear, ordinary differential equations with respect to y 
with constant coefficients 

[[/]){*} = { o} (i9) 


with the boundary conditions 

( 20 ) 

There is one set of Eqs. (19) for region A, shown in Fig. 
2, and another set for region B. Thus, solutions are 
obtained for each region as described below, and the 
solutions are connected by the inter-region continuity 
conditions that are obtained from Eq. (20). (The inter- 
region continuity conditions are discussed in more detail 
subsequently). The solution to the homogeneous set of 
Eqs. (19) is obtained by letting 

{ X } = ( 21 ) 


where, for convenience, the following transformations 
are introduced for the width coordinate y 


y = y - b for y > b 

y = b - y for y < b 


( 22 ) 


The vertex at y = b (Fig. 2b) is then defined by y = 0 , 
in both regions, and y is always positive. Substitution 
of Eq. (21) into the homogeneous differential equations 
leads to the eigenvalue problem 

[VV]A. 4 + [VH 2 + [U] = 0 (23) 

which can be written in the standard form 

[D] {<p} = !{<p} (24) 


where [D] is an unsymmetric matrix given by 


[Z)] = [C 0 ]-'[C 1 ] 


~[U] 

[0] 

-1 

IV] 

[W] 

_[0] 

[/]. 


-[/] 

[0] 


(25) 


The matrices [C Q ] and [Cj] have dimension 6n + 1 
by 6n + 1 , and the matrices [ U] and [ V] have 
dimension An + 1 by 4n + 1 . The matrix [ W] , which 
is associated with the coefficient matrix [W] of the 
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operator d 4 /dy 4 , has dimension 2 n by 4n + I and has 
nonzero elements only in the top 2 n by 2 n submatrix 
(in rows and columns corresponding to the functions F k 
and G k ). The zero matrices [0] have dimension In 
by 4/i H- 1 and the identity matrices [/] have dimen- 
sion In by 2 n . 

The eigenvalue problem has been solved by the 
highly efficient combination shift QZ algorithm. 12 For a 
general laminate, there is a total of 6n + 1 eigenvalues 
A, 2 ; thus there are 6n + 1 pairs of eigenvalues X , that 
are either real or complex. One component of the pair 
has a positive real part and the other component of the 
pair has a negative real part. For the case of a cross-ply 
laminate, the elastic compliance coefficients p 26 , P 36 
and P 45 and S j6 are zero in all layers, and the system of 
equations (Eq. (19)) decouples into two systems of equa- 
tions for the stress functions O and 4* . The system as- 
sociated with O determines the solution for the in-plane 
problem (plane strain), and the system associated with 
'F determines the solution for a twisting load. The so- 
lution spaces for the eigenvalue problems associated 
with these two systems have dimension 8 n and 4n + 2 , 
respectively. Eigenvectors associated with the eigenval- 
ues are determined from the homogeneous form of Eq. 
(19). The total solution is obtained as a superposition of 
the eigenvectors multiplied by generalized coordinates 
(complementary solution) plus a particular solution. A 
particular solution to Eq. (19) can be obtained by assum- 
ing {X (y) } is constant. A particular solution was ob- 
tained from the far-field classical lamination theory 
stresses, and is outlined in the Appendix. 

As mentioned previously, the stiffened panel prob- 
lem is solved by dividing the structure into two regions. 
If independent stress functions are assumed for the two 
regions, two uncoupled eigenvalue problems are ob- 
tained when the principle of minimum complementary 
energy is applied. There are n interfaces in region A 
and, for a general laminate, 6n + 1 pairs of eigenvalues 
contribute to the solution. There are n s interfaces in re- 
gion B and hence 6^ + 1 pairs of eigenvalues. Express- 
ing the solutions in these regions in terms of the 
eigenvectors and the particular solutions gives 

6 n + 1 

{*“} = £ a f {tp.}e X ^+ {X°} 

<“> 

{A = X bj{Xj}e***+ { Xf } 

j= i 

where the superscripts a and b denote the regions, and 
the subscript p denotes the particular solution. In the 
expansions in Eq. (26), only negative eigenvalues are 
considered since local coordinates y have been estab- 
lished with their origins at the thickness discontinuity 
vertex, and stresses decay with distance from this point. 
The local coordinates were chosen in this manner to 


avoid potential numerical overflow problems associated 
with the evaluation of positive exponential functions. 
The constants a . and b. are determined by imposing 
the free edge conditions along the flange edge AB and at 
the vertex B, and the inter-region continuity conditions 
along BC (see Fig. 2). The free-surface conditions 
require 

{ Y a (y = 0) } = 0 k = 1 -> n f (27) 


where = n- n s * The continuity conditions between 
the regions are 


{ Y a ( j = 0) } = {^(>’ = 0)} 


m -^3 + [ 01-^2 + [*]-£+ [ S ] 

, dy s dy 2 dy 


( 


+ [5] 

, dy 3 dy 2 dy 


O'" 0 = 0)} 


{K b (3> = 0)} 


k = 


1 +/ 


Results 


(28) 


The approximate analytical solution presented 
above has been implemented in a FORTRAN code for 
general laminates. Input to the local analysis program 
is simple and concise, and includes loading conditions, a 
description of the geometry, stacking sequence informa- 
tion, and ply material properties. Plies or layers may 
have different thicknesses and they may be of different 
materials. The loading conditions for the local analysis 
are far-field in-plane stresses that are determined from 
classical lamination theory based on force and moment 
resultants obtained from a global analysis. Output from 
the code includes the eigenvalues, eigenfunctions, and 
constants in the eigenfunction expansion. Once the 
eigenvalues, the eigenfunctions and the constants in the 
eigenfunction expansion are determined, transverse and 
through-the-thickness stress distributions can be easily 
obtained at any specified location. 

Results are presented for two examples. The accu- 
racy of the present analysis of the first example is as- 
sessed by comparing the present analysis results with 
results from an elasticity solution 2,3 and results from a fi- 
nite element analysis. Computer execution times are 
also presented for this example. Further, a brief discus- 
sion of the influence of the modeling detail used in the 
present analysis on the interlaminar stress predictions is 
presented. The second example illustrates the breadth of 
information that can be obtained from the present analy- 
sis. Stress contour plots, and transverse and through-the- 
thickness stress distributions near the thickness disconti- 
nuity are presented. The interaction between the 
through-the-thickness and the transverse stress gradients 
in a region near the thickness discontinuity is also dis- 
cussed. 
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Orthotropic Skin and Flange Plates 

The first example application of the analysis is a 
stiffened graphite-epoxy plate subjected to pure bending 
by a moment M y applied about the x-axis. The stiffener 
flange is a ' rectangular, 8-ply quasi-isotropic 
[±45/0/90] laminated plate as indicated in the upper 
sketch of Fig. 4. The skin is an 8-ply orthotropic 
[±45/90 2 ] laminated plate that is also rectangular, but 
is wider than the flange. Each ply in both the flange and 
skin plates has a uniform thickness f (*) = 0.005 in. and 
moduli in the planes of elastic symmetry given by 

E n = 20 msi E 22 = £33 = 2.1 msi 

Gp = Gp = G 2 3 = 0.85 msi (29) 


where 1 , 2, and 3 refer to fiber, transverse, and thickness 
directions, respectively. To consistently compare results 
of the present solution with results of an elasticity solu- 
tion, 2 ’ 3 each of the laminates is replaced by an equivalent 
homogeneous orthotropic plate with smeared compli- 
ances obtained from the respective [A] extensional 
stiffness matrices (see Table 1). The present analysis, 
however, is not subjected to this limitation; that is, it is 
possible to model each layer separately. Furthermore, 
the layers may be monoclinic, and are not limited to spe- 
cially orthotropic layers. 

Typical results for the interlaminar shear stress a 
and the interlaminar normal stress C zz along the skin- 
stiffener interface obtained from the present analysis, 
from an elasticity solution and from a finite element anal- 
ysis of the stiffened plate are shown in Figs. 4a and 4b, re- 
spectively. The interlaminar shear stress a for this 
problem is zero valued, as a result of the uncoupling of 
the in-plane problem from the torsion problem for spe- 
cially orthotropic material systems. The elasticity solu- 
tion is a tedious solution based on complex-variable 
respresentations of the stress functions and a boundary 
collocation method. 2 The finite element analysis was 
conducted using the displacement based finite element 
code CLFE2D. 13 The code CLFE2D was developed for 
analyzing two-dimensional generalized plane deforma- 


tion boundary value problems. Thus, only a cross-sec- 
tion of the stiffened plate was modeled in the finite 
element analysis. The finite element model of the thick- 
ness discontinuity region used in the present study has 
4033 nodes distributed over 1.5 inches of length. Re- 
sults from the analytical model presented herein were 
obtained using 16 layers to model each of the skin and 
flange plates, and yields a total of 1 87 eigenvalues in the 
thick region (region A), and 91 eigenvalues in the thin 
region (region B). The stresses shown in Figs. 4a and 4b 
are normalized by M y t % /21 and are plotted as a function 
of the normalized distance y = y/t from the flange 
termination vertex, where M y , t s and 1 denote the ap- 
plied moment, the skin plate thickness, and the moment 
of inertia of the skin plate, respectively. The solid lines 
correspond to the present analysis, the dashed lines cor- 
respond to the elasticity analysis, and the symbols corre- 
spond to the finite element analysis. 

The finite element predictions shown in Fig. 4 are 
discontinuous at the flange-skin interface, and cannot be 
obtained exactly at the interface without interpolation. 
Therefore, to provide the most accurate stress data, finite 
element predictions are plotted for Gauss point locations 
in the flange and skin plates, just above and below the in- 
terface, respectively. Only one value is given for the in- 
terface data computed by the elasticity solution and the 
present analysis, since these analyses satisfy stress con- 
tinuity across the interface. As the figures indicate, the 
interlaminar normal stress predictions of all three meth- 
ods agree very well. The shear stress predictions of the 
current method and the elasticity solution also agree 
very well, but there is some discrepancy in the finite el- 
ement predictions near the vertex of the thickness dis- 
continuity. The correlation between the finite element 
results, and the present analysis results and the elasticity 
solution is expected to improve as the finite element 
mesh is refined. In addition, the results indicate that a 
significant boundary layer is present in the region near 
the thickness discontinuity. 

Cohen and Hyer 2,3 have indicated that it is possible 
for elastic singularities to exist at the thickness disconti- 
nuity vertex. The current analysis does not include a sin- 
gularity explicitly in the formulation. However, singular 


TABLE 1. Smeared Elastic Properties for Flange and Skin Plates 


Laminate 

E, 

(msi) 

Ey 

(msi) 

E z 

(msi) 

G XZ 

(msi) 

G yz 

(msi) 

G xy 

(msi) 

^XZ 

(msi) 

V yz 

(msi) 

v *y 

(msi) 

Flange [±45/0/90] s 

8.05 

8.05 

2.16 

0.85 

0.85 

3.09 

0.164 

0.164 

0.303 

Skin [±45/90 2 ] 

3.83 

11.71 

3.83 

0.85 

0.85 

3.09 

0.176 

0.107 

0.201 
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behavior is predicted by the present analysis. The singu- 
lar behavior becomes apparent for the example problem 
by increasing the number of mathematical subdivisions 
used to model the flange and skin plates. The interlami- 
nar shear and normal stress predictions along the skin- 
stiffener interface are shown in Figs. 5a and 5b, respec- 
tively, for solutions where 4, 8 and 16 equal thickness su- 
blayers are used to model each of the flange and skin 
plates. The solid lines correspond to the 4-layer model, 
the dashed lines correspond to the 8-layer model, the 
dashed-dotted lines correspond to the 1 6-layer model, and 
the circular symbols correspond to an elasticity solution. 
The present analysis results based on the 16-layer model 
agree more closely with the exact elasticity solution in the 
very local neighborhood of the flange termination at 
y = 0 than with the corresponding results based on 8- 
layer and 4-layer models. This observation is particularly 
true for the interlaminar shear stress o results. The in- 
terlaminar shear stress is equal to its prescribed value of 
zero along y = 0 for all of the discretization levels 
shown in Fig. 5. Very steep gradients are encountered in 
the transverse interlaminar shear stress distribution as the 
flange termination vertex is approached. As the number 
of layers in the analytical model increases, the predicted 
maximum stress increases, and the location of this maxi- 
mum stress shifts toward the vertex at y = 0 , thereby ap- 
proaching the elasticity solution. The interlaminar 
normal stress, in contrast, is not prescribed along y = 0 
but is singular at the vertex, and the predicted stress at this 
point increases monotonically as the number of mathe- 
matical layers in the model increases. However, the 
present analysis results away from the vertex are nearly 
identical for all of the discretization levels that are shown. 
A similar observation is made for the interlaminar shear 
stress. These observations are significant from a design 
perspective, in that depending on the design and failure 
criteria used in conjunction with the stress analysis, the 
results based on the 4-layer model may be as adequate for 
predicting failure as the results based on the 16-layer 


model, and the 4-layer model takes much less time to ex- 
ecute. Execution times recorded on a Personal Iris 4D / 
35 Workstation for models based on five through-the- 
thickness discretizations are provided in Table 2. These 
execution times indicate that doubling the number of lay- 
ers in the analytical model increases the run time by an 
order of magnitude. 

Another method for improving solution accuracy, 
while keeping computer execution times as low as possi- 
ble, is to adjust the relative thicknesses of the sublayers 
used to model the plates rather than increasing the num- 
ber of sublayers used. 14 This approach is basically 
equivalent to mesh refinement. Thick layers can be used 
for the regions away from the skin-flange interface, and 
an increased number of thin layers is used near the inter- 
face where the gradients are the most pronounced. Inter- 
laminar shear stresses along the skin-stiffener interface 
are shown in Fig 6. These results were obtained from a 
model in which the flange and skin plates are represented 
by two mathematical sublayers. Stress predictions are 
presented in Fig. 6 as a function of the ratio ty /t 2 , where 
ty is the thickness of the mathematical layer used to 
model a portion of the flange or skin plate adjacent to the 
skin-stiffener interface, and t 2 is the thickness of the 
mathematical layer used to model the remaining thick- 
ness of the plate. As the thickness of the layer adjacent 
to the skin-stiffener interface is reduced, the results for 
O yz from the present analysis approach the elasticity so- 
lution. Similar trends have been found for the corre- 
sponding interlaminar normal stress o zz results. 

Isotropic Pl ates Subjected to Uniform Extension 

The second example is the analysis of a stiffened 
aluminum skin panel subjected to an extensional load 
Ny applied to the skin panel. The configuration of the 
stiffened skin is the same as for the previous example. 
Both the flange and the skin are 0.04 inches thick and 
have a modulus of elasticity E = 10xl0 6 psi and a 
Poisson’s ratio v = 0.30 . There is no material -property 
mismatch between the flange and skin plates. Therefore, 


TABLE 2. Execution Times for Five Through-the-Thickness Discretizations 


Number of 
layers 

Total number 
of eigenvalues 

Number of skin 
plate 

eigenvalues 

/(*) 

(inches) 

Time 

(cpu seconds) 

4 

43 

19 

0.010 

2.21 

5 

55 

25 

0.008 

3.82 

8 

91 

43 

0.005 

11.7 

10 

115 

55 

0.004 

26.0 

16 

187 

91 

0.002 

100 
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this example shows the effect of the geometric disconti- 
nuity at the flange termination on the stress field, and 
provides a simple case for studying the interactions be- 
tween the transverse and through-the-thickness gradients 
in the stress field near the thickness discontinuity. 

Distributions of the interlaminar shear stress a 
and the interlaminar normal stress a^ 7 obtained from the 
present analysis are shown in Fig. 7. The stress compo- 
nents a and a xz are zero for this loading condition 
and material system. Stresses are normalized by the 
skin-panel far- field stress a . There is a stress concen- 
tration at the flange termination point, where the stresses 
have a lobe-shaped distribution near the corner of the 
flange. Large through-the-thickness gradients and in- 
plane gradients in both G^ and G 77 are observed near 
the flange termination point (y = 0 ). Away from the 
flange termination point, the interlaminar stresses are 
equal to zero. 

The through-the-thickness and in-plane gradients 
of the stress components, and their interactions, are 
shown explicitly by plots of the stress distributions along 
the lines z/t s = 0.0, z/t s = -0.025 and 
z/t s = -0.098 shown in Fig. 8 by curves A-A, B-B and 
C-C, respectively, and at y/t s = 0.022 (section D-D), 
shown in Fig. 9. When the flange termination region is 
approached in the direction y , gradients in G yy develop, 
as shown in Fig. 8, as load is transferred from the stiff- 
ened to the unstiffened plate regions. These transverse 
gradients in o yy are equilibrated by steep through-the- 
thickness gradients in a , as shown in Fig. 9. The 
through-the-thickness o yz distribution is self-equilibrat- 
ing; that is, it integrates to zero through the thickness (its 
resultant force is zero valued). A through-the-thickness 
distribution of the interlaminar normal stress G zz is also 
shown in Fig. 9. The interlaminar normal stress o zz 
changes sign abruptly as the flange-skin interface is ap- 
proached; it changes from compression in the flange to 
tension in the skin. There is also a local variation in the 
stress distribution as the thickness discontinuity is ap- 
proached. This local variation is consistent with the con- 
tour plots in Fig. 7 and is caused by the different rates of 
change of a with changes in y at different locations 
through the thickness of the comer region. 

Concluding Remarks 

A geometrically linear approximate analysis meth- 
od for predicting complete three-dimensional stress 
fields in laminated composite structures with local thick- 
ness discontinuities has been presented. The analysis 
method uses stresses obtained from a global two-dimen- 
sional analysis as traction boundary conditions for a local 
approximate three-dimensional analysis to determine 
stresses in regions with local thickness discontinuities. 
The analysis assumes that a state of generalized plane de- 
formation exists. The approximate local analysis is 
based on the principle of minimum complementary ener- 


gy and uses stress functions with assumed through-the- 
thickness functional dependence. Input to the local anal- 
ysis is simple and concise, and includes loading condi- 
tions obtained from the global analysis, a description of 
the geometry, stacking sequence information, and ply 
material properties. Interfacial and through-the-thick- 
ness stress distributions are easily obtained at any speci- 
fied location, and are not limited to the stiffener-flange 
interface or interfaces between plies. 

The present analysis method has been applied here- 
in to determine the stresses in the flange termination re- 
gions of a stiffened laminated composite plate and a 
stiffened aluminum plate subjected to simple bending 
and transverse extension loading conditions, respective- 
ly. Stress predictions from the present analysis have 
been compared with corresponding finite element results 
and results from an elasticity solution, and shown to 
agree very well. The correlation is shown to improve as 
the number of layers used in the analytical model is in- 
creased. Furthermore, it has been demonstrated that, al- 
though the analysis model does not explicitly include a 
singularity, a singularity appears in the predictions as in- 
creasing stresses at the flange-skin interface with in- 
creasing number of layers in the analytical model. In 
addition, changing the relative thicknesses of the mathe- 
matical layers used to model the skin and flange ele- 
ments of a stiffened plate was shown to provide a 
modeling discretization approach for improving stress 
predictions without increasing the number of layers and 
the computational requirements. The computational ef- 
ficiency and accuracy of the method make it well suited 
for use in preliminary design and optimization studies of 
advanced composite structures with local structural de- 
tails. 
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A ppendix 

Elements of particular solution vectors 


6 2 *(» 0 0 
-6 -4AO) 0 0 

-6 4A( 2 > 6 2 A< 2 > 
6 -2/i< 2 > -6 -4A< 2 > 


a 2 </> (/«(») 2 


a 2 < 2 >(M 2 )) 2 

m ] ( a <2) ) 2 



( 30 ) 


where the subscripts t and b denote the top and bottom 
of a layer, respectively, the tilde denotes classical lami- 
nation theory stresses, and n f is the number of layers in 
the flange plate. The particular solutions for L k and H k 
in the flange plate are obtained by solving the system of 
equations: 







V 

-2 

-*<*> 

0 

0 . . 



0 

1 

0 

0 . . 


i 

L, 

2 

0 

-2 

-h (2) . . 


2 

0 

0 

0 

1 . . 


h 2 






H . 




d£>ft< 2 > 


(31) 


[ dl n b f ) h (n f )] 

The particular solutions in the skin plate are given in 
Eqs. (32) and (33). 


The particular solutions for F k and G k in the 
flange plate are obtained by solving the system of equa- 
tions: 
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Figure 1. Panel with thickness discontinuity. 



(a) Applied loads (b) Cross-section model 

Figure 2. Cross-section model and applied loads. 



Figure 3. Generic layer. 
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skin 




(a) Interlaminar shear stress 



(b) Interlaminar normal stress 


Figure 4. Comparison of interlaminar stress distributions along the skin-stiffener interface predicted by 
the present analysis, an elasticity solution, and a finite element analysis. 
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(b) Interlaminar normal stress 

Figure 5. Interlaminar stress distributions calculated from an elasticity solution 
and the present analysis using three through-the-thickness refinements. 
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Figure 6. Influence of layer discretization on the interlaminar shear stress 
predictions along the skin-stiffener interface. 



(a) Interlaminar shear stress U (b) Interlaminar normal stress G , 

yz zz 

Figure 7. Interlaminar shear and normal stress distributions near a flange termination point of 
stiffened aluminum panel subjected to an in-plane skin load N y 


0.16 

- 0.12 
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Figure 8. Distribution of transverse in-plane stress at three through-the-thickness 
locations of a stiffened aluminum panel subjected to an in-plane skin load N y . 


t, 






Figure 9. Through-the-thickness interlaminar shear and normal stress distributions at section 
(D-D) of a stiffened aluminum panel subjected to an in-plane skin load Ny, 
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